Lasing at half the Josephson frequency with exponentially long coherence times 
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We describe a superconducting device capable of producing laser light in the visible range at half 
the Josephson generation frequency, with the optical phase of the light locked to the superconducting 
phase difference. An earlier proposed device, the so called 'Half- Josephson Laser (HJL) [Phys. Rev. 
Lett. 107, 73901 (2011)], cannot provide long coherence times, because of spontaneous switchings 
between the emitter states. To circumvent this we consider N 3> 1 emitters driving an optical 
resonator mode. We derive a general model that captures essential physics of such devices while 
not depending on specific microscopic details. We find the conditions under which the coherence 
times are exponentially long, thus surpassing the fundamental limitation on the coherence times of 
common lasers. For this we study the noise in the device. In particular, we are interested in the 
rate of large fluctuations of the light field in the limit where the typical fluctuations are small. The 
large fluctuations are responsible for switching of the laser between stable states of radiation and 
therefore determine the coherence time. 

PACS numbers: 42.55.Px, 85.25.Cp, 74.45. +c, 42.60.Mi 



All numerous degrees of freedom of a bulk supercon- 
ductor are frozen at low temperature. So its state is 
characterized by a single variable, the quantum phase of 
the superconducting order parameter—. In principle this 
phase remains constant, that is coherent, for infinitely 
long time. Another example of a highly coherent system 
is a laser. There a macroscopic number of photons form 
a coherent state characterized by an optical phase and 
an amplitude 2 . In contrast to superconductors the co- 
herence time is not infinite, but limited by the intrinsic 
noise in the device, caused by spontaneous photon emis- 
sion. This leads to a phase drift and thus to coherence 
loss at typical timescales of2> Td oc — n/T, with n the num- 
ber of photons in the mode and T the cavity escape rate. 
Coupling these two very different systems may open up 
novel ways of controlling either of them, allowing for ma- 
nipulation of one with the other. 

Recently, a device has been proposed^ that com- 
bines the coherence of laser light with the coherence of 
the superconducting order parameter, dubbed the 'Half- 
Josephson Laser' (HJL). The light in this device is gener- 
ated by a biased semiconductor p-n junction realized on 
a nanowire. To acquire a narrow emission spectrum of 
the junction, quantum dots are defined at the p-n inter- 
face. The coupling of light to the superconducting phase 
is achieved by connecting both ends of the nanowire to 
superconducting leads, thus inducing proximity gaps at 
the quantum dot regions. Even though the device is bi- 
ased with a voltage much larger than the superconduct- 
ing gap, any direct charge transfer is effectively blocked 
by the potential barriers in the p-n junction. Instead, 
electrons and holes jump to the quantum dots and there 
recombine and emit a photon. This might happen in two 
different ways. A "blue" photon at about the Joseph- 
son frequency, ujj, is emitted when Cooper pairs from 
either side of the junction recombine. The phase of this 
"blue" photon is locked to the superconducting phase 



difference. A process resembling spontaneous paramet- 
ric down-conversion^ allows for the emission of two "red" 
photons in stead of a single "blue" one, at about half the 
Josephson frequency, uij/2. The setup as described up to 
here was proposed earlier as the 'Josephson Light Emit- 
ting Diode' 6 (JoLED). It was shown that it can be used 
for quantum manipulation purposes^. The JoLED is also 
the basis for the setup of the HJL. The latter exploits the 
emission of "red" photons by the JoLED and amplifies it 
by stimulated emission in an optical resonator. In con- 
trast to the case of the JoLED, where "red" photons are 
emitted spontaneously, the "red" photons in the reso- 
nant mode of the HJL form a coherent state. The phase 
of this coherent state is locked to the superconducting 
phase difference. 

Lasing was found for the "red" emission with the opti- 
cal phase locked to the superconducting phase difference, 
where only two values of the optical phase are allowed, 
having a phase difference of tt. The importance of the 
phase locking in this device is that it can result in expo- 
nentially long laser coherence times. The idea to prolong 
coherence times by phase locking has been thoroughly 
explored in the field of lasers^. Using, for instance, the 
technique of injection locking^ one could pump a laser 
with another one with longer coherence time, so that the 
pumped laser inherits this longer coherence time. How- 
ever, in that case the coherence time is still limited by the 
fundamental noise processes in the pump laser, with some 
coherence time Td cc - The essential difference with the 
HJL is that the laser light inherits the coherence time of 
the superconductors, which is infinitely long. It will not 
become infinitely long, though, because the fundamental 
noise in the laser is still present and causes switchings 
to the lasing state with opposite phase, with exponen- 
tially small probability. Hence the coherence times can 
be exponentially long. 

Unfortunately, despite this promising prospect, the 
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HJL does not exhibit these exponentially long coherence 
times because switchings between quantum dot states re- 
sult in a shorter coherence time. The magnitude and 
phase of the laser light depends on the actual quantum 
state of the JoLED in the HJL. Not every state supports 
lasing however, and the ones that do give rise to vari- 
ous magnitudes and phases of the laser light. Switchings 
between the quantum dot states therefore cause deco- 
herence of the HJL. The switching rates are suppressed 
because they occur in combination with the emission of 
off-resonant photons, which is allowed because of the fi- 
nite linewidth of the optical resonator. For an optimal 
choice of parameters it is found that, despite this very dif- 
fer decoherence mechanism, the coherence time is, Td cc , 
similar to that of the common laser and much shorter 
than the envisioned exponentially long coherence times. 

The design for the HJL is based on a single quantum 
emitter, putting rather high demands on the quality fac- 
tor of the optical resonator. It was proposed^ to relax 
these demands by considering a large number, N ^> 1, 
emitters in a single resonator, increasing proportionally 
the upper bound on the damping rate at which lasing 
is achieved. As an added benefit, intensity fluctuations 
are expected to reduce as 1/yN, thus realizing a less 
fluctuating light intensity. 

In this article, we thoroughly investigate the alterna- 
tive idea, where, instead of a single JoLED, we consider 
a large number of emitters forming a dipole moment to 
drive the optical resonator mode. The model that we 
introduce will be quite general as we will only consider 
macroscopic quantities and leave out the microscopic de- 
tails. The crucial aspect of the design, the coupling of 
the quantum states to both superconducting leads, will 
of course remain, since it both drives the laser and in- 
duces the phase-lock between optical phase and super- 
conducting phase difference. The high demands on the 
quality factor of the resonator can be relaxed because it 
is driven by a big number of emitters. 

We find that stationary lasing is possible with two 
equivalent values of the optical phase, as in the HJL, 
and with exponentially long coherence times owing to 
the phase-lock. Rather surprisingly, the sole increase in 
number of emitters does not directly guarantee the in- 
creased number of photons in the mode. The average 
dipole moment of all emitters will turn out to be zero 
when the emitters are coupled to the resonator mode 
only. Incoherent emission to other modes is required to 
create a population imbalance^ resulting in a nonzero 
average dipole moment. Owing to the large number of 
emitters, single switchings will not cause decoherence as 
in the HJL. Rather, they induce minor fluctuations of 
the dipole moment, transferring to detuning fluctuations 
in the resonator mode. These again transfer to fluctu- 
ations of magnitude and phase of the lasing state. The 
detuning fluctuations compete with the intrinsic fluctua- 
tions of the resonator mode. When both are sufficiently 
small, the superconducting phase difference will remain 
locked to a single value of the optical phase. The major 



decoherence mechanism of the laser is then provided by 
large fluctuations, represented by the tails of the noise 
distributions. These bring the optical field to the state 
with opposite phase. Large switchings occur rarely, on 
exponentially long timescales. The importance of such 
a feature is hard to underestimate. Indeed, examining 
these processes we have found coherence times (surpass- 
ing the fundamental limit for common lasers), scaling as 
Tdcc ~ exp(n), up to some critical value n s ^> 1, where 
it saturates. While the linewidth of the laser is directly 
related to this coherence time, there are large tails in the 
line shape associated to the noise in the dipole moment 
as well as to the intrinsic quantum noise. 

In the model we study, we neglect the noise source 
causing decoherence of the superconducting phase differ- 
ence: voltage fluctuations. The superconducting phase 
difference is fixed only for ideal voltage bias. In case 
of nonidcal bias, the voltage fluctuations will drive the 
phase evolution and cause the drift of the superconduct- 
ing phase difference. The corresponding coherence time 
can be estimated using the voltage correlation function 
for a junction with impedance Z 

(V(t)V(t')) = 6(t - t')k B TZ (1) 

with ksT the thermal energy for temperature T. With 
this the average value of the phase is given by 

/ e i<t>(t)\ = e -^«0(t)0(*')» = e -±£k B TZt , 2 s 

where the second Josephson relation and the voltage cor- 
relation function were used for the final expression. From 
this we find a coherence time of Td cc ,sc = h 2 /4e 2 kBTZ. 
The coherence time can thus be prolonged by reducing 
temperature or impedance of the nanostructure. Alter- 
natively one can increase the number of Josephson junc- 
tions to reduce the noise. The latter approach is used 
in metrology to realize the Josephson voltage standard. 
A voltage drift of a few nV per hour can be realized-^ 
in that way, leading to coherence times on the order of 
10 12 seconds. Hence the voltage fluctuations may be dis- 
regarded when studying the decoherence in the device. 

This article is organized as follows. In Sec. U we derive 
the governing equations for a broader class of devices, 
characterized by a large number of emitters. We start 
by briefly treating the model of the HJL, as described in 
Ref. 0, and subsequently argue that the resulting equa- 
tions can be generalized to represent a wide range of de- 
vices. In Sec. [II] it is shown that lasing with a phase- lock 
occurs for these general equations in certain parameter 
regimes. There are distinguished regimes, separated by 
first and second order phase transitions. In Sec. [TTT] con- 
ditions are derived to determine the parameter regime 
where the noise can be considered in linear approxima- 
tion. The frequency dependence of the noise, when in a 
stationary lasing state, is calculated. The noise in lin- 
ear approximation cannot cause the decoherence of the 
laser owing to the phase lock. However, occasional large 
fluctuations will cause decoherence. These are treated in 
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Sec. lIVl for which the formalism of optimal paths is used. 
Finally, we conclude in Sec. [V] 

I. MODEL 

The purpose of this section is to derive a rather gen- 
eral but simple model describing a novel class of lasing 
devices, those are driven by superconducting leads bi- 
ased at voltage V , have the optical resonance frequency 
at half the Josephson generation frequency, 2eV/H, and 
contain many quantum states capable of emitting light. 
The purpose of the general approach is to capture the 
essential properties of this class of devices, without in- 
volving specific microscopic details of the used model. In 
this section we start by briefly repeating the essentials of 
the model for the HJL at microscopic level and extend 
it to contain a large number of effective emitters, yet 
at the microscopic level. We formulate the generic fea- 
tures of the model which are not specific to the details 
of the microscopic description. Finally, we present the 
general but simple model where all these generic features 
are incorporated. We formulate the model in terms of a 
Fokker-Planck equation for three variables. 

A. The HJL with a single emitter 

The HJL is described in Ref. as an optical resonator 
driven by a single quantum emitter. The emitter is a 
double quantum dot in a p-n semiconductor nanowire, 
capable of containing up to two electrons and two heavy 
holes. The energy of the corresponding discrete quantum 
states is determined by orbital and Coulomb energies and 
is represented by Hamiltonian H QD . Owing to the su- 
perconducting proximity effect, a gap is induced at the 
quantum dots for both electrons and holes, represented 
by Hamiltonian 

ffso = A*c t CL + Ahhfhi + H.c. (3) 

where the A e are the proximity induced pair poten- 
tials, with the phases 4> e .h retaining the values of the 
corresponding superconducting leads. Owing to gauge 
invariance, only the superconducting phase difference, 
= 4>e — 4>h will have a physical significance. Interac- 
tion of the double quantum dot with the resonator mode 
occurs according to 

H pb = hub% + G(tfx + bx f ) (4) 

where u> = ujq — eV/h is the frequency detuning from res- 
onance, with u>o the resonance frequency of the optical 
resonator and V the bias over the device. The operator 
x = + ft-fCj,, taken in rotating frame, represents the 
electron hole recombination required for the photon emis- 
sion. It can be seen as the dipole moment operator that 
drives the resonator mode. Coefficient G is the coupling 
strength of the dipole moment to the resonator mode. It 



is given by the dipole strength owing to quantum fluc- 
tuations of the electric field in the resonator mode. The 
interaction Hamiltonian is assumed to be perturbativc 
as to allow for a large photon number before nonlincari- 
ties become manifest, which will induce saturation of the 
dipole moment. 

The dynamics of the resonator mode can be treated in 
a semiclassical fashion, because a large number of pho- 
tons will occupy the mode. In view of this, the b\ b and 
b'b will be replaced by their respective expectation val- 
ues in the photon-dependent part of the full Hamiltonian. 
A scaled quantity is defined, A = G(b). Within this ap- 
proximation, the expectation value of the dipole moment, 
x m (X) = (m\x\m), is found by taking the derivative of 
the spectrum E m (X) to A*, with E m (X) and eigenstates 
|m) associated to the full Hamiltonian, H QD + H sc + H ph . 
Hence, the value of the dipole moment, x m , depends on 
the particular eigenstate of the double quantum dot. On 
the one hand, the dipole moment depends on the optical 
field, A, but on the other hand, the evolution of A is de- 
termined by the dipole moment. This leads to a set of 
self-consistent equations 

A = - liw + - I A - i— x m {X), x m = -QT^-. (5) 

Stationary lasing is found when A = for finite values of 
A. This state must be stable against small perturbations. 

To get lasing, the rate at which photons are created 
in the resonator mode, 2(G 2 /S)Im[x m (A)/A], needs to 
exceed the cavity escape rate, T. In common lasers, the 
imaginary part of the susceptibility, Im[x m (A)/A], results 
from pumping to excited states. Usually, to achieve las- 
ing, the pumping must be sufficiently intense to invert the 
population. For the HJL, it is the superconductivity that 
creates the required complex susceptibility so that the 
population invertion is not necessary. In fact, the ability 
of the superconductors to absorb or supply Cooper pairs 
undermines the notions of groundstate and excited state 
in the quantum dots. For instance, the induced pair po- 
tentials enable the possibility of transitions from lower to 
higher energy states, while emitting a photon, with the 
required energy provided by the bias. This feature allows 
for cycles of emission or absorption of photons, in order 
to (de)populate the resonator mode. 

Decoherence in the HJL with a single emitter is caused 
by spontaneous switchings between quantum dot eigen- 
states. Owing to the finite cavity lincwidth, emission 
of off-resonant photons occurs inducing these switchings 
They occur at a slow rate, r sw ~ T/n <C T, and therefore 
introduce extra dynamics in the HJL described by a mas- 
ter equation. This is in addition to the already present 
dynamics of Eq. It was found in Ref. that even 
for an optimal choice of parameters, not every quantum 
dot eigenstate couples to the resonator mode strongly 
enough to support lasing. The ones that do support las- 
ing will yield different stationary values A s . Since only 
some states support lasing, the optical field may even ex- 
tinguish after a switch, to turn on again after another 
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one. Hence, the switchings change X s causing large in- 
tensity and phase fluctuations. Most notably, for every 
lasing state X s there is another one at — A s . When going 
from some lasing state to a nonlasing state and return- 
ing to the same lasing state again, the optical field can 
go to both ±A S with equal probability. Hence the switch- 
ings are responsible for decoherence. The coherence time 
is determined by the switching rate, Td oc — n/T, and is 
similar to that of the common laser—. 

B. The HJL with many emitters 

Even though the HJL, as described above, has a rather 
specific design, the framework provided by Eqs. ©-(lU 
is quite general. It applies to a class of devices where a 
biased Josephson junction contains a structure of which 
the eigenstates couple to both superconducting leads and 
which can emit light by electron hole recombination. The 
light is emitted into an optical resonator. When the 
structure in the Josephson junction is large enough it 
can be divided into N 3> 1 blocks. Every block repre- 
sents an effective quantum emitter, which is a quantum 
system with a number of discrete eigenstates that can 
emit light. Every emitter must couple to both super- 
conducting leads and must emit light into the resonator 
mode. In principle the emitter interact with each other, 
but we may model them to be independent. As an exam- 
ple of such a structure, one might consider an ensemble 
of p-n semiconductor nanowires with a double quantum 
dot, the nanowire that was used in the previous subsec- 
tion, or one might consider a two dimensional electron 
gas between the superconducting electrodes^. 

The dipole moment corresponding to the structure in 
the Josephson junction is determined by its energy, which 
is taken in a rotating frame of reference. When the struc- 
ture contains N emitters the total energy is given by the 
sum of the energies of the emitters. Assuming that emit- 
ter i is in a state with energy Ei, the total energy is sim- 
ply E — J^i Ei- This yields a value of the total dipole 
moment given by D = dE/d\*, which should be inserted 
in Eq. ((SJ) to look for solutions that support lasing, X s . 

Every emitter is subjected to spontaneous switchings 
that change its eigenstate. Therefore, at timescales 
longer than the switching times it is not the total en- 
ergy (dipole moment) that is relevant but the average 
total energy (dipole moment). The switching dynamics 
defines a master equation for each of the emitters, where 
the probability p l m of emitter i to be in state m evolves 
as 

= E i W U S nm ~ KM* (6) 

with W^-n the transition rate from state \n)i — > \m)i 
in emitter i. The quantity in brackets [• • • ] defines a 
transition matrix, which necessarily has a null vector 
to ensure the conservation of probability. The null vec- 
tor gives the stationary occupation probabilities, (p l m ) s . 



With these, the average total energy is given by (E) = 
Y,i, m (Pln)sE l m , where E l m the energy of state m in emit- 
ter i. Because of the coupling to the resonator mode, the 
energy depends on A and the total average dipole mo- 
ment is (D) = d(E)/dX*. Now this quantity determines 
the lasing solutions of Eq. ([5]) . 



C. The average value of the dipole moment 

Let us discuss a complication that might easily arise 
when considering this model or designing a device. In 
our approach we consider spontaneous switchings only 
to be a result of off-resonant photon emission. In the 
limit of a large number of emitters the associated transi- 
tion rate matrix, Eq. (j6|), turns out to be approximately 
symmetric, with W^. n ~ W^. m . The transition rates 
Wjjj.jj are calculated using Fermi's Golden Rule. They 
are, up to some constants, given by the product of a 
matrix clement squared, |(TO(A)|_ff ph |n(A))| 2 , for a given 
emitter i, and a Lorentzian shaped density of states, 
piElJ = (l/27r)hT /{(hT) 2 + (E* mn n with E) nn the 
energy difference between the eigenstates \m(X))i and 
\n(X))i. We note that the photon number generically re- 
mains the same during a transition; A is constant. Only 
after the transition may relaxation to the new station- 
ary state occur. For a large number of emitters however, 
a single switching is not enough to change the electrical 
field in the mode significantly. The stationary state A s re- 
mains almost the same after the transition and therefore 
also the eigenstates and energies remain approximately 
the same. Consequently, the transition rates also remain 
the same. Since the energies E\ nn have a large mismatch 
with the resonance frequency (the emitted photons are 
far off-resonant) they do not depend significantly on de- 
tuning, cj, and we have E l mn = E l nm HT, for fixed 
A. The transition rate from state |m(A))j to \n(X))i is 
therefore equal to the transition rate of the reversed pro- 
cess. Hence, because the transition rates remain approx- 
imately the same after a transition, the transition rate 
matrix is approximately symmetric. 

Innocent as it may seem, this feature of the transition 
matrix poses a problem for our device, because as a re- 
sult the average value of the dipole moment, (D), must 
be approximately zero. The null vector for a symmetric 
transition matrix must be a vector with all entries equal. 
Therefore, in the stationary solution to the master equa- 
tions, all the states in an effective emitter i have the same 
probability p\ to be occupied. The total average dipole 
moment is then given by (D) = J2iPlJ2m x m- But for 
the summation over m we find 

m m 

where Tr^ is the trace over the space of states of emitter 
i. The equality to zero follows because Tii[H] is indepen- 
dent of A*. Consequently (D) — 0. We note that this is 
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generic for any kind of system governed by our Hamilto- 
nian and where the relevant states defining the dipoles, 
are equally populated. Naively one would expect that a 
large number of emitters would result in a large number 
of photons in the resonator mode. We see that this is not 
so. On the contrary, the total dipole moment could be 
zero and no photons are present. 

We need to adjust the model used in Ref. to get a 
nonzero average dipole moment. This is done by making 
the transition rate matrix asymmetric. Then the occupa- 
tion probabilities for the states of an emitter will not be 
equal anymore, resulting in an imbalance in populations 
of states. There is no a priori reason to favor, for instance, 
the state with higher energy over the state with lower 
energy, for higher occupation probability. The notion of 
higher or lower energy of states has become obscured by 
the superconducting pair potentials. Hence, the created 
imbalance is not the same as a population invertion. Im- 
portantly, recombination of electron-hole pairs will still 
be characterized by the phase-lock, because they are still 
coherent with the superconductors. Despite introducing 
incoherent transition rates the special characteristics of 
the device are thus maintained. 

An asymmetry in the transition rate matrix can come 
about in various ways. The simplest way is to take into 
account the photon emission into the environment, rather 
than to the resonant mode. If the rate of this incoherent 
process is comparable with the emission rate to the reso- 
nant mode, the average dipole moment is of the order of 
its maximum value per emitter. In a realistic system we 
expect relaxation transitions between the emitter states 
that do not involve photon emission. These relaxation 
channels will generally yield asymmetric transition rate 
matrices. One can also engineer transition rates to get 
a sufficiently asymmetric transition matrix. A straight- 
forward approach is to change the rates that are already 
present by changing the electromagnetic environment. A 
switching to a state with higher (lower) energy is ac- 
companied by emission of an off-resonant photon, with 
energy eV - E nm (eV + E nm ), where E nm > T. By 
changing the density of states in favor of one of these 
two, the spontaneous emission rate for the correspond- 
ing process increases with respect to the other one, thus 
creating the desired imbalance in transition rates. This 
might be done by inserting the device in a low quality 
factor resonator with resonance frequency about cither 
of the frequencies. The low quality factor will ensure 
the rapid escape of the photon, before it is reabsorbed 
again. Alternatively, the electromagnetic environment 
can be changed by shining with a red or blue detuned 
laser on the device, favoring transitions to lower and 
higher energy eigenstates respectively, proportional to 
the laser intensity. This resembles the technique of re- 
solved sideband cooling/heatingH. Another possibility is 
to add new transition rates to the system by adding non- 
superconducting leads to the device and directly pump 
electrons and/or holes into the optically active region. 
This however leads to big challenges in the technical de- 
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D. Phase dependence of the energy 

We have discussed the general features of possible mi- 
croscopic realizations of a HJL with a large number of 
emitters. Let us construct a model that does not rely on 
microscopic details. The relevant quantity to deal with 
is the A-dependent part of the total energy of the device, 
which determines the dipole moment. We show that the 
energy has a generic dependence on the phase difference 
between the superconducting leads. 

Generally, we note that both common lasers and super- 
conductors are systems characterized by a spontaneously 
broken {/(l)-symmetry, a phase symmetry. As such they 
are characterized by a order parameter which has a mag- 
nitude and a phase. Any actual realization of the order 
parameter bears no apparent connection to the original 
unbroken symmetry. However, as a remnant of the sym- 
metry, every value of the phase is equally probable and 
corresponding to states with equal energy. States with 
different phases are degenerate so that the energy of the 
system is independent of the phase. 

Our model combines both these systems, lasers and su- 
perconductors, which leads to a phase lock between the 
optical phase and superconducting phase difference. In- 
deed, owing to the interaction between these two systems 
a phase dependence must enter the energy of the system 
so that not both of the spontaneously broken symme- 
tries can be maintained. The emission of a single pho- 
ton is caused by the recombination of an electron-hole 
pair, with the latter carrying the phase of their respec- 
tive condensates. The phase of the terms in the relevant 
interaction Hamiltonian, Eq. Q, is therefore given by 
4>\ — <M/2, with 4>x the optical phase and 4>a the super- 
conducting phase difference. In fact, the full interaction 
depends on the cosine of this phase combination, lead- 
ing to a similar dependence in the energy. As a conse- 
quence an energy minimum is found for a specific value 
of <p\ — 0a/2, which is a manifestation of the phase-lock 
between optical phase and superconducting phase differ- 
ence. A further inspection of the phase dependence of 
the energy shows that there must in fact be two equiva- 
lent values of the phases. Changing <j)& — > a + 2-7T leaves 
the system invariant since the superconducting phase dif- 
ference has not changed, but the phase in the interac- 
tion has changed to <j>\ — ir — 0a /2, so we conclude that 
E(cf)\) = E((j)\ — 7r). Hence, the original spontaneously 
broken [/(l)-symmetries of the Hamiltonian have been 
reduced to a single one and a spontaneously broken sym- 
metry under rotations of 7r of the optical phase. The 
remaining J7(l)-symmetry comprises invariance under a 
simultaneous change in optical phase and superconduct- 
ing phase difference such that the quantity 4>\ — a /2 is 
unchanged. We stress that this reasoning applies to any 
system that combines two order parameters with sponta- 
neously broken [/(l)-symmetries, with an interaction as 
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in Eq. Q. Any such system manifests a phase-lock. 
E. Energy 

To find lasing in the general model, we need to have 
an expression for the energy of the device. This allows 
us to extract essential physics from the model. 

Generally, the energy should depend on the value of 
light field in the resonator, b = (6), and the pair poten- 
tials, A e , A/j, in the superconductors: E(b, A e , A h ). We 
will use parameter limits where the four variables are rel- 
atively small, enabling us to use a Taylor expansion for 
the energy. This limit corresponds to a weak drive of 
the resonator mode. To define 'small' we go back to the 
model where the dipole medium consists of many emit- 
ters. We can associate an energy scale, E s , to each of 
these emitters. The induced pair potentials should be 
small compared to this scale, A e; h <C E s . In a realistic 
device this can be engineered by making suitable (tunnel) 
barriers between the superconductors and emitters. For 
the second limit, it is reasonable to take the interaction 
between dipole and resonator mode in the weak coupling 
limit. Hence the coupling constant G of the interaction 
Hamiltonian, Eq. 21 is much smaller than E s . Therefore, 
this interaction perturbs the dipole energy by a value of 
the order of E s only when a large number of photons, 
«o = \bo\ 2 , occupies the resonator mode. We want to 
make sure that the stationary value of photon number, 
n s , is always much smaller than no, so that the use of an 
expansion of the energy in terms of b is always justified. 

When writing down the Taylor expansion of the energy, 
E(b, A e , A/j), only a limited number of terms are allowed 
because of the specific phase dependence that is required. 
As was discussed in the previous subsection, the energy 
should depend on the phase combination, ip — 2<j)\ — ((f) e — 
<f>h). Therefore the lowest order phase dependent term 
in the energy must be proportional to A* A/jfe 2 +complex 
conjugate. This is a fourth order term in the Taylor 
expansion. All other terms of equal or lower order are 
proportional to powers of |6| 2 or |A e ,/j| 2 . Writing down 
only the terms with immediate relevance, the energy is, 
up to fourth order in Taylor expansion 

E = E + Q'\b\ 2 + V|fo| 4 + ^{b* 2 + b 2 ) (8) 

In this expression the phase of the induced pair potentials 
was absorbed into the optical phase, so it is redefined as 
4>\ — > ip/2 — <fi\ — 0a/2. Furthermore, we absorbed the 
induced pair potentials in the prefactor A = C|A e A^|, 
with C a constant. The primes in Q' and Q" refer re- 
spectively to first and second derivatives at \b\ 2 = of 
the function 51(|6| 2 ), which is independent of the induced 
pair potentials. The lowest order term, Eq, is a sum of 
f2(0) and some other terms depending on | A e; /j | 2 but not 
on b and b* . These terms will therefore not be relevant for 
the dipole moment. There are also terms proportional to 
|A e ,fc| 2 |&| 2 , but these are small compared to O', because 



they are of higher order in the Taylor expansion, and 
can safely be neglected. It is necessary to keep the \b\ 4 
term in the expansion, since this nonlinear term stabilizes 
a possible lasing instability. Without this term the light 
field would grow exponentially under unstable conditions. 
When calculating the dipole moment later on, we will see 
that fi' adds up to the detuning from resonance. Hence, 
the value of fl' can implicitly be changed by changing 
the detuning. In fact, to get lasing the detuning should 
be such that ft' effectively becomes small, of the order 
of A or smaller, or goes to zero. Because the detuning 
M < eV/K, this also sets the scale \Cl'\ < eV/%. We 
stress that this is an effective change of Q' . In the expan- 
sion of Eq. © the SI' term is still considered to be larger 
than the fourth order terms. With effectively 57' < A, 
the \b\ 2 term is of the same order as, or smaller than the 
|6| 4 term so the latter should be maintained in the ex- 
pression for the energy. With the above expression for 
the energy we have managed to separate the dependence 
of the energy on b and A e ,/j from the microscopic details 
of the driving medium. The latter are all contained in 
the prefactors of the Taylor expansion. 

Substituting Eq. © into Eq. ([3]), we find an explicit, 
general self-consistency equation. Noting that A = G(b) 
and dE/db* = G(d), we have 

b= - yiu+7£j b-i(Ab* + n"\b\ 2 b), (9) 

where h = 1 and Co = uj + f2' is the redefined detuning. 
This equation describes the deterministic evolution of the 
electric field, 6, in the resonator mode. We note that an 
equivalent complex conjugate equation exists. 



F. Fokker-Planck equation 

The energy was considered to be a fixed quantity in 
the above section, corresponding to the average value of 
the dipole moment. To fully account for all dynamics 
in the HJL with many emitters, we also have to consider 
noise. This can be incorporated in a single Fokker-Planck 
equation that encompasses two noise sources. The first 
noise source is the traditional quantum noise due to pho- 
ton emission from the mode 1 ^. The second noise source 
describes the fluctuations of parameters in Eq. (|8|) . 

To understand the latter let us note that eq (|8]) is valid 
not only in average, but also for each specific configu- 
ration of the emitter quantum states. The parameters 
D,', il" and A are configuration dependent. Since the 
emitters undergo spontaneous transitions the parameters 
fluctuate accordingly. Explicitly taking into account all 
the details of the fluctuations would greatly complicate 
the model. Instead, we employ the central limit theo- 
rem stating that any ensemble of equal distributions can 
be described as a Gaussian distribution when the num- 
ber of distributions goes to infinity. Therefore, since we 
have a large number of emitters, as described above, we 



7 



may regard the parameters as stochastical quantities with 
Gaussian distributions. 

In principle the noise in all parameters is expected to 
be of the same order of magnitude relative to their aver- 
age values. But since f2' is of lower order in the Taylor 
expansion of Eq. (jSJ , its fluctuations are more important 
than that of Vt" and A. We will therefore only consider 
a single noise source by taking VL" and A to be equal to 
their average values, while Q! d' e + (l, with fl' e being 
the average value and f2 being the deviation. The latter 
appears in Eq. Q as an addition to the detuning. The 
noise can therefore be seen as fluctuations in detuning or 
in the resonant frequency. The corresponding variance, 
(i7 2 ), is an important parameter of our general model. 
We shall take into account that emitters switch on a cer- 
tain timescale that determines autocorrelations of these 
detuning fluctuations. The autocorrelation function is 
thus given byA5, 



«n(t)fi(i + r))) - (Q 2 )e-Tl T l, 



(10) 



where 7 is the inverse autocorrelation time of the emitter 
ensemble 

It is possible to capture all relevant dynamics com- 
pactly in a single nonlinear Fokker-Planck equation. This 
is an equation for the probability distribution function 
P(b, Cl, t) of three real values, two of them are repre- 
sented in the optical field and the third in the detuning 
deviation. The drift and diffusion terms in the equation 
respectively describe deterministic and noise driven, in- 
deterministic processes^. The Fokker-Planck equation 
for our device encompasses two Gaussian noise sources, 
describing fluctuations in three quantities. First, there is 
the dynamics of b which is subjected to a noise source, b, 
describing the fundamental fluctuations in b. This is the 
quantum noise associated to the discrete changes of pho- 
ton number in the resonator mode. It satisfies, (6(f)) = 
and (b(t)b*(t')) = TS(t — t'). Second, as discussed, the 
noise in Q! is equivalent to detuning fluctuations. This is 
represented by the quantity J7, with average (Q) = and 
variance (Cl 2 ). The Fokker-Planck equation reads 



dP(b,n,t) 
dt 



= 2Re 



_5 

db 



i(u + n"\b\ 2 + n) + - 



bP + iAb*P 



r d 2 P 
2 dbdb* 



B - - d 2 P 

JLnp + (n 2 )^J- 
an on 2 



(ii) 



The single derivative terms on the right-hand side of 
Eq. (fTTj) are the drift terms. Evaluating (b) — J dbbP by 
using Eq. (fTTj) . will yield the self-consistency equation (|9]) 
for (b). We remark that in this context the brackets, 
(•••), denote the average over the classical distribution 
P, and not the expectation value of the operator b. The 
diffusion is described by the second derivative terms on 
the right-hand side of Eq. (|TT|) . 



Let us note that all the parameters we use have the 
dimension of frequency. It is convenient to make them 
dimensionless measuring them in relative units of T/2. 
This amounts to setting r — > 2 in Eq. (|lip . For clarity, we 
rescale the detuning deviation in Eq. pip to appropriate 
units, by taking x = Q/y/A 2 — 1, which is the scale of 
the frequency domain where lasing is possible, as we will 
show in Sec.Ql] Additionally we have (x 2 ) = (fl 2 )/(A 2 — 
1). As a consequence, the term in Eq. ([lip proportional 



to iilbP will change to i\J A 2 
their form. 



1 xbP- Other terms retain 



II. LASING 

In the previous section we have derived a model for 
a class of lasers that is driven by superconductivity and 
that exhibits a phase-lock. Here we will show that lasing 
with a phase-lock is indeed possible and we give expres- 
sions for the lasing threshold and the number of photons 
in the resonator mode. We find three different regimes, 
with a different number of stationary lasing states. 

A full description of the dynamics of lasing would re- 
quire solving Fokker-Planck equation (jTTJ) , but this is 
difficult, if not impossible, even for a relatively simple 
system such as ours. Instead, we will first search for sta- 
tionary lasing solutions to the deterministic equation of 
motion, Eq. ©, thus disregarding noise. Once these so- 
lutions are found we will study the effects of noise when 
near a stable stationary lasing state. This is the topic of 
the subsequent two sections. 

To find the desired stationary solutions to Eq. ©, with 
6 = 0, we rewrite them as real equations. Defining b = 
x + iy leads to a set of equations 



We conclude this section by stressing that this Fokker- 
Planck equation applies to a whole class of lasers, where 
a biased Josephson junction contains a structure of which 
the eigenstates couple to both superconducting leads and 
which can emit light by electron hole recombination into 
the resonator mode. 



x = (u>-A + Sl"n)y, y = -(u> + A + tt"n)x (12) 

where n — x 2 + y 2 is the photon number in the resonator 
mode. These equations can be readily solved to find an 
expression for n and for the phase of 6, which was defined 
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FIG. 1. Lasing in the device. The upper plot shows the num- 
ber of photons, n, for the stationary solutions of Eq. (|12|l . 
as a function of the rescaled detuning \, with Q." < 0. Full 
(dashed) lines represent stable (unstable) solutions, labeled 
by n±. Three parameter regimes, (I, II and III), are found 
based on the number of stationary solutions. The vertical 
dotted lines represent boundaries between these regimes. Be- 
low this plot is a figure that illustrates the existing stationary 
solutions in the x-y plane, for each regime. Solutions with 
n > occur twofold with opposite sign. Filled (empty) cir- 
cles represent stable (unstable) stationary states. The lower 
figure shows the dependence of the regimes on the parameters 
\A\ and X- The diagonal lines are the lasing thresholds \± 
as discussed in the main text. The line between regimes II 
and III is the one where n± becomes real-valued. The plot of 
the upper figure is taken along the dashed line. When cross- 
ing the boundary along arrow(s) (i) and (ii) [(iii)] , the system 
undergoes a second [first] order phase transition. 



as tp/2 = (f>x - </>a/2, 



-(f) 



xM 2 



Q' 1 



[±i - X] , 



(13) 



where we conveniently rescaled the detuning, u — 
xV A 2 — 1. Hence, we find stationary lasing solutions and 
a phase-lock. Note that the equations are invariant under 



the transformation x — > —x and y — > — y, implying the 
occurrence of two solutions which differ in phase by 7r, as 
predicted in the previous section. Indeed, Eqs. (fTB"]) are 
invariant under the transformation ip±/2 — » <p±/2 + ir. 
Hence n± and tp± represent four stationary solutions. 
Besides these, there is another stationary solution, being 
x = y = 0, delivering a total of five possible stationary 
solutions. We remark that tp± only depends on A and 
not on x an d Furthermore, a large photon number 
can be achieved when Q," <C A. 

Studying the stationary solutions of Eq. (fT3"|) we find 
that three different regimes that differ in the number of 
physical solutions. Since n± represents the number of 
photons in the resonator mode, it must be either a real 
and positive quantity or it must be zero. As a result 
we note that lasing solutions, with n± > 0, can only be 
achieved when A 2 > 1. Then, depending on x, three 
different regimes exist (Fig. Q]). We choose J7" < but 
similar regimes can be found when VL" > 0. For x < and 
|x| > 1, both n± < so that the only physical solution 
is n — 0. We call the parameter regime with a single 
solution regime I. Regime II is given by the condition 
that |x| < 1, where n+ < and n_ > 0. Including 
n = there are now three physical solutions. The third 
regime is when x > and \x\ > 1, where both n± > 0, 
corresponding to a total of five physical solutions. 

Even though there can be multiple physical solutions, 
not all of them will support stationary lasing, because 
some will be unstable against small perturbations. To 
find the condition for stability we expand Eq. 0, in 
terms of x and y, up to first order about a stationary 
value. This yields 



d_ f8x 
dt Uy 



■l + 2Cl"xy 



| + 2n"y 2 
-1 - 2n"xy 



(14) 



where x and y are at the stationary values, 5x and Sy are 
the respective deviations from the stationary values and 
Eq. (|12j) was used to rewrite some terms. When an eigen- 
value of the matrix is negative (positive), a perturbation 
in the direction of the corresponding eigenvector decays 
(grows) exponentially to (from) the stationary value of 
x and y. The case where both eigenvalues are negative 
corresponds to a stable stationary lasing state. After a 
fluctuation in any direction, the system will relax to the 
stationary state again. With one positive and one nega- 
tive eigenvalue there is a stable and an unstable direction 
against perturbations, thus defining a saddle point. The 
case where both eigenvalues are positive cannot occur, 
because the trace of the matrix is smaller than zero. The 
condition of the stability therefore depends on the sign of 
the determinant of the matrix, positive being the stable 
stationary state. Evaluating the determinant and again 
using Eq. (fTiZj) we derive that the stability for a physical 
solution n± > is determined by sign [4 (A 2 — 1)(1 =p %)] , 
whereas the stability of the solution at n = is deter- 
mined by sign [(A 2 — 1)(1 — x 2 )] ■ 

To illustrate how the system relaxes to a stable sta- 
tionary lasing state after a small perturbation, we give 
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the eigenvalues of the matrix in Eq. (fT4]) . These are, for 
the case of n_, 771,2 = — 1 ± i\J^{A 2 — 1)(1 + x) — 1 and, 
for the case of n +1 rji^ = — 1 ± i^jA(A 2 — 1)(1 — x) — 1. 
Hence, when A is sufficiently large, the system oscillates 
about the stationary point with a frequency ± ~ 1/A 
and relaxes to it in ~ A/2ir rotations. 

For f2" < the stability of the stationary solutions in 
the various regimes is as depicted in figure [T] Regime 
II is the most relevant one for lasing as the n — so- 
lution is unstable whereas the n_-solution is stable. In 
regime III the n_- and n = 0-solutions are all stable, 
whereas the n+-solutions are unstable. When Q," > 0, 
the qualitative picture is the same. Quantitatively, the 
location of the regimes changes. In the phase diagram 
of Fig. [TJ the regimes I and III are interchanged so that 
the plots for f2" > are mirror images of the plots for 
£1" < 0. Additionally, the n + (n_) solution is now the 
stable (unstable) one. 

The boundaries between the regimes can be seen as las- 
ing thresholds because they separate the lasing regimes 
from the nonlasing regime. The boundaries between 
regime I and II and between regime II and III are lines 
in parameter space where the determinant of the ma- 
trix of Eq. (|14[) is zero at n = 0. They are given by 
X± = ±1. For negative O" the value x- (x+) corre- 
sponds to the boundary between regime I and II (II and 
III). The boundary between regime I and III is deter- 
mined by the condition A = 1. At this boundary, the 
n = solution is still stable, but two pairs of stable and 
unstable solutions arise at values n > 0. When gradu- 
ally crossing the boundaries from I to II or II to III, the 
system undergoes a second order phase transition [ar- 
rows (i) and (ii) in Fig. Q]. Thereby the n — 0-state 
becomes unstable. In contrast to this, when crossing the 
boundary from regime I to III [arrow (iii) in Fig. [T] the 
system will undergo a first order phase transition, be- 
cause all physical stationary solutions will arise directly 
at nonzero values, two stable and two unstable. We note 
that A cannot be increased indefinitely as the pairing po- 
tentials must remain small enough to satisfy the adopted 
approximations. 

To conclude this section, we find that regime II is the 
only regime where lasing occurs without an interruption 
caused by a large fluctuation. This is because only there 
the 6 = solution is unstable. In contrast to this, in 
regime III a large fluctuation can bring the device from a 
lasing state to the stable nonlasing state at n = thereby 
interrupting the emission of light. Regime II occurs in 
a small interval at large detuning, where \J A 2 — l|x| = 
|w| = \uj + Cl'\ ~ \A\ < |ST|. The number of photons 
in the resonator mode, being of the order of |-A/f2"|, is 
large when |0"| <C \A\. In the remainder of this work, 
we will assume we have chosen the parameters A and x 
such that the system is in regime II. We take Vt" < so 
that the n_ solution will be the lasing solution. 
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FIG. 2. Noise and the noise spectrum in the HJL with many 
emitters. On the left-hand side is an illustration of fluctu- 
ations in the optical field, b — x + iy, where we neglect for 
clarity the detuning noise on the z-axis. The value of b evolves 
along the lines. The dashed line, indicated by A, represents a 
large fluctuation, where 6, starting in the vicinity of a stable 
lasing point, crosses the unstable point to the other stable 
lasing point. This process shall occur on exponentially long 
timescales, Ti arg(! (see main text). The solid erratic lines, in- 
dicated with B, represent small fluctuations in the vicinity 
of the stable points. These fluctuations typically occur on a 
timescale r sma ii (see main text). The two kinds of fluctuation 
contribute independently to the noise spectrum, of which an 
illustration is shown at the right-hand side. The sharp peak, 
indicated with A, corresponds to the large fluctuations. It 
is a Lorentzian shaped peak with exponentially small width. 
The small fluctuations provide a broad-band background, in- 
dicated by B. Details on the peaks can be found in the main 
text. 



III. NOISE AND ITS FREQUENCY 
DEPENDENCE 

Let us consider fluctuations in the device. It is im- 
portant to note that because of the phase-lock the laser 
cannot lose coherence by a random drift of the optical 
phase as in common lasers. On the contrary, after a 
small fluctuation a damping force returns the system to 
the stable stationary lasing state. The average value of 
the optical phase is therefore retained. As in the HJL 4 , 
the only possibility to lose coherence is when the device 
switches to another stable stationary lasing state. This 
requires a sufficiently large fluctuation in photon number 
or optical phase. We envision these large fluctuations to 
occur only rarely, on exponentially long timescales, im- 
plying that the noise in b must be linear, i.e. dominated 
by small fluctuations about the stationary lasing state. 
In this section we investigate the HJL in the linear noise 
approximation. There are two noise sources causing in- 
trinsic quantum fluctuations in the resonator mode and 
fluctuations in detuning. Therefore, the steady state of 
the HJL can be described by a multivariate Gaussian 
probability distribution, for which we calculate the vari- 
ances. Several parameter limits will be discussed con- 
cisely for the photon number variance. We also find an 
expression for the corresponding noise spectrum. The 
large fluctuations that cause decoherence of the laser will 
be the topic of the next section. 

We first describe an intuitive picture of the shape of 
the full noise spectrum, containing the spectrum of both 
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large and small fluctuations. A sketch of this is given in 
Fig. [5] The large fluctuations occur at exponentially long 
timescales, T largc , satisfying ln(r largl ./r amall ) ^> 1. Here, as 
we will see later in this section, the small fluctuations oc- 
cur at much shorter timescales, determined by the largest 
of the response time of the resonator mode, ~ 1 /V, or the 
autocorrelation time of the detuning noise, ~ 1/7, yield- 
ing r amall = max(l/r, I/7). Because of this separation 
of timescales, both kinds of fluctuation contribute inde- 
pendently to the noise spectrum. This is given by the 
Fourier transform of the correlator of b(t), 



S(iz) 



d(t -t') ((b*(t)b(t')))e 



iv(t-t') 



(15) 



At timescales r largo the correlator is expected to decay 
because of the switching to another lasing state. The 
switching changes the optical phase but retains the sta- 
tionary photon number n s = \b s \ 2 . The correlator decays 
as 



((b*(t)b(t>)))=n s 



|*-t'|/n a 



The noise spectrum of this correlator is a Lorentzian with 
exponentially small width, 1/T largc and with peak height 
£(0) = 2n s T largo . At timescales r small small fluctuations 
have a dominating contribution to the correlator. Writ- 
ing b(t) = b s + a(t) for a small fluctuation, a(t) <C b s , 
we have ((b*(t)b(t'))) = (a*(t)a(f)). The requirement 
that the fluctuations a(t) are small implies that they are 
Gaussian distributed as will be shown below. One ex- 
pects the noise spectrum to consist of three peaks. One 
peak is associated with the fluctuations in detuning, and 
has the width ~ 7. Two side peaks, shifted with respect 
to zero detuning, are associated with the intrinsic quan- 
tum fluctuations and have width ~ 1. These peaks are 
shifted because of the oscillatory behaviour of small devi- 
ations from equilibrium position, which was described in 
the previous section in the contex t of stability analysis . 
In fact, the shifts should be ±^/4(A 2 - 1)(1 + x) - 1. 
The noise spectrum approximately corresponds to a sum 
of three Lorentzians. Here the peak heights must be 
5(0) ~ 2(|a(0)| 2 )r Bmall < 2n s r largo . Adding up the two 
contributions, the total noise spectrum is a high and very 
narrow peak on a background formed by three low and 
broad peaks, as depicted in Fig.O In this section we will 
concentrate on the small fluctuations that are responsi- 
ble for the broad background in the noise spectrum of 
the figure. 

A full description of the small fluctuations can be found 
by linearizing the Fokker-Planck equation of Eq. (TTTT) . 
This amounts to expanding b to b(t) — b s + a(t) and 
only retaining up to linear terms in a, a* and \. The 
resulting Fokker-Planck equation can be solved exactly^, 
but we will not do this, because we are only interested 
in the steady state probability distribution. The latter 
is fully characterized by the variances of the variables of 
the system. 

Instead of directly using the linearized Fokker-Planck 
equation of Eq. (fTTj) to calculate the variances, we use the 



equivalent system of Langevin equations for a set of ap- 
propriate variables. The advantage is that in the course 
of this procedure, we directly obtain the noise spectrum 
associated to the linear noise. The linearized form of 
Eq. (|11[) is fully equivalent^ to the following system of 
Langevin equations, 



a = - i{^J A 2 - lx + 2n"n s ) + 1 



-i(A + n"b 2 s ) a* - ix(t)VA 2 - 1 b s + &(t), (16) 

X = -7* + &(*)■ 

Here it is noted that the first equation is complex and rep- 
resents two equations. There are two different Langevin 
forces satisfying 

(&(*)> =0, (&*(W)> = <*(*-*'), (17) 

with other correlators being zero. The source £&(£) rep- 
resents the quantum noise in the resonator mode, while 
the source (^(t)} represents the noise in the detuning. 
The current set of equations is not the most convenient 
one because the variables a and a* are not orthogonal to 
each other. It will be most natural to work in a system 
with a variable Sn, representing radial displacements or 
photon number changes, and a variable Sip, represent- 
ing azimuthal displacements or phase changes. The re- 
lation between these two bases comes about in the fol- 
lowing way. One can write the deviation a = Sx + iSy 
in vector form in a Cartesian coordinate system as a = 
Sx x + Sy y. Alternatively, in a cylindrical coordinate 
system a — Sr f + ^/ri^Sip ip, where the ^fn~ s follows be- 
cause a is taken about b s , which is separated from the 
origin with a distance ^/nl. The variations Sip repre- 
sent the variations in optical phase, but the variations 
Sr correspond to the fluctuations of the square root of n. 
Hence, Sr — Sn/(2y/n^). Using the transformation rela- 
tions, f — cos tp x + sin ip y and tp = — sin tp x + cos ip y, 
one can derive expressions for Sx and Sy to find that 



Sn 



+ iy/n^Sip 



Using this we can rewrite the above set of Langevin equa- 
tions. Reminding the reader that we use for the lasing 
solution n s = n_ and tp s = tp- [Eq. (1131) ]. we arrive at 



Sn = -4n s V(^ 2 -l) Sip + £„(t) 

5<p = -25ip + — A 2 - 1(1 + x )Sn (18) 



Here we have three Langevin forces satisfying 
(&,(*)) =0, (Ut)Ut'))=n s S(t-t'), 
(ip{t)) = (UW)) = ~ t'), (19) 

<&(«)> =0 (^(t)^(t'))=2j(x 2 )6(t-t>), 
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The other correlators are zero. 

With the change of variables, we now need to find 
the steady state probability distribution to the Fokker- 
Planck equation for the variables Sn, Sip and \. This is 
given byi£ 



P s (6r,S(p,x) 



1 



1 —T <— — 1 — 



(20) 



V(27r) 3 Det(S) 

(6n,5<p,x) 

((Sn 2 ) (SnSip) (Snx)\ 

(SnS^) (Sip 2 ) (5<px) 

(Snx) (Sipx) (X 2 ) J 



where (x) is simply the variance of the distribution 
of the detuning fluctuations. Now using the Langevin 
equations we can find the variance matrix S according 
to the following procedure: first, we take the Fourier 
transform of Eqs. (fT5|) . For any variable g(t) we have 
g(t) = J du <7„exp[— if]. The resulting system of equa- 
tions is solved for 5n v , Sip v and x„. To get the variances, 
we calculate the inverse Fourier transform at t = of the 
expectation value of products of the Fourier transformed 
variables. For instance, we have 



(SnSip) 



1 

2^ 



dv(Sn v Sip v ) 



(21) 



According to this procedure, we find the variances and 
write them in a form that will be convenient later on 



(Sn 2 



2( 7 + 2)(A 2 -l)(x 2 ) 



(l+x)[4(A 2 -l)(l + X )+7 2 + 2 7 ] 



+ 



2n s 



A 2 



vV> 

(SnScp) 
n s 

(fax) 

n s 

(<fy>x) 



(A 2 -l)(l + x ) 
l(A 2 -l)(x 2 ) 



2[4(yl 2 -l)(l + X )+7 2 + 27] 
1 



8n s 



An s ^(A 2 -l) 

4(A 2 -l)(x 2 ) 
4(A 2 -l)(l + x)+7 2 + 27 

- 7 VI 2 ^T(x 2 ) 
4(A 2 -l)(l+x)+7 2 + 27 



(22) 



We remark that some variances diverge when approach- 
ing the lasing threshold, which is when \A\ — > 1 or 
X —> — 1. Also, the combination 4(^4 2 — 1)(1 +x) 1S equal 
to the determinant of the stability matrix of Eq. (fT4|) . A 
larger value of the determinant implies higher stability, 
as the variances decrease. 

With the expressions for the variances, we are able to 
quantify in which parameter regimes the linear Fokker- 
Planck equation and the steady state distribution of 
Eq. (|20|) are valid. Their validity is ensured when Sn <C 
n s and 8<p <C 1, so that all expressions in Eq. (|22[) must 
be much smaller than one. Note that under these con- 
ditions, there is no a priori bound on x which is already 



purely linear. It will be particularly useful to study the 
expression for (Sn 2 ) /n 2 in a bit more detail, because it is 
always of the order of, or larger than the other ones. Ad- 
ditionally, this is an interesting quantity, because large 
fluctuations in Sn lead to switchings and are thus impor- 
tant for the long coherence times. 

The variance of the fluctuations in photon number is 
studied in a few limits. For reasons of simplicity and 
clarity, it will be most convenient to be sufficiently far 
from the lasing threshold, so we take x "C 1 and A> 1. 
Keeping this in mind, the first limit to be studied is the 
one where A 3> 7 . Then we have 
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n 2 2 17 



2)(X 2 
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(23) 



We note a few properties. To start with, the term pro- 
portional to (x 2 ) is associated to the noise in the detun- 
ing, whereas the l/n s term is associated to the intrinsic 
quantum noise of the resonator mode. The latter exactly 
corresponds to the fundamental noise expected in any 
common laser 2 . In line with this division, there are two 
distinct regimes when considering the dependence on n s . 
For small enough values of n s the intrinsic quantum noise 
dominates so that the variance decreases inversely pro- 
portional to n s . Increasing n s beyond a critical value the 
variance saturates as the l/n s term becomes unimpor- 
tant and the noise in detuning starts to dominate. We 
further note that the magnitude of the fluctuations is in- 
dependent of A and independent of 7 in the limit 7 <C 1 . 
The criterion to have small fluctuations now leads to the 
requirements that n s 3> 1 and (7 + 2)(x 2 ) <C 1. 

The two previously mentioned limits of 7 in the de- 
tuning noise dominated regime can be intuitively under- 
stood. If 7 <C 1, it is smaller than the decay rate of 
the resonator mode so that the fluctuations in x are slow 
enough for the optical field to adiabatically follow them. 
The number fluctuations are then purely determined by 
the equilibrium dipole distribution. When 7 > 1, the 
optical field cannot follow the dipole fluctuations adia- 
batically. However because A ^> 7, the optical field will 
respond fast to the dipole fluctuation before it fades away. 
In fact, the optical field will start to rotate about the 
stationary point with a period t ~ 2ir/A and the fluc- 
tuation only fades away after ~ A/2n rotations about 
the stationary lasing state. The variance of the num- 
ber fluctuations is now determined by the instantaneous, 
out-of-equilibrium value of the fluctuation of x, as deter- 
mined by £%(t) in Eq. (TIT)) , leading to an extra factor of 
7 compared to the previous case with 7 <C 1. 

We continue by studying the variance of photon num- 
ber in the second limit, which is A <C 7 and leads to 
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(24) 



In this limit, again there is in the dependence on n s a 
division in two regimes: one where the variance is domi- 
nated by the noise in detuning and one where the intrin- 
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sic quantum noise of the resonator mode is most impor- 
tant. In contrast to previous case there is now a depen- 
dence on both A and 7. Here, the criterion to have small 
fluctuations leads to the requirements that n s 3> 1 and 

A 2 {x 2 ) < 7- 

There is also an intuitive way to understand the detun- 
ing noise dominated regime in the limit ^4^7. Here, 
the dynamics of the detuning fluctuations is on much 
shorter timescales, I/7, than the dynamics of the res- 



[{v + VD~T)2 + - v^D^T) 2 + 1] 



The term proportional to (x 2 ) is directly related to the 
noise in the detuning, whereas the other term is directly 
related to the quantum noise of the resonator mode. 
Which of the two terms is the most important, depends 
of course on the parameter regime as was discuss above. 
The spectrum is not separable into a sum of Lorentzians, 
but can be approximated by it, particularly when D ^> 1. 
The noise spectrum has a central peak at v$ = and 
there are two sideband peaks at v± = ±\/D — 1. As- 
suming again that A 1 and X <C 1, the central peak 
resembles a Lorentzian peak with width 27 when 7 <C A. 
For 7 > 4, the central peak almost vanishes compared 
to the sideband peaks. The latter are well described 
by Lorentzians with width 2. These peaks indeed cor- 
respond to the intuitive picture discussed above. 

We have now described the linear fluctuations in the 
HJL and derived conditions on parameters for which fluc- 
tuations are small, (Sn 2 )/n 2 <C 1. The noise in various 
parameter regimes was studied. For small n s , but still 
much larger than 1, the intrinsic quantum noise of the 
resonator mode was found to have a dominating con- 
tribution to the photon number fluctuations, while for 
large n s , the noise in the detuning is the most important. 
There, the magnitude of number fluctuations depends on 
7 and the ratio A/7. The noise spectrum corresponding 
to the small fluctuations was found and has a central 
peak with varying width, and two sideband peaks with 
width 2. 



IV. LARGE FLUCTUATIONS 

In the previous section we have studied small fluctua- 
tions in the HJL under the assumption that they domi- 
nate the fluctuation spectrum. These small fluctuations 
do not cause decoherence of the laser. Decoherence is 
caused by large fluctuations that switch the stationary 
lasing state to one with a phase difference of ir, see fig[T] 
In this section we will study the large fluctuations thor- 
oughly to find the timescales at which the switchings 



onator mode, 1/A, so that the latter has a reduced time 
window to respond to the detuning fluctuations. Com- 
pared to the earlier case with 1 <C 7 <C A, the reduction 
factor is 4A 2 /~/ 2 . 

As part of the calculation done to find the variances, 
we can also find an expression for the noise spectrum, 
which is defined as S(v) = {\5a u \ 2 ) = (\5n„\ 2 ) / j (An s ) + 
n s {\5tp v \ 2 ). Defining D = A(A 2 - 1)(1 + x) we have 



(25) 

I 

would occur. This gives the decoherence time of the HJL. 

The probability to have a large fluctuation can be cal- 
culated using the formalism of optimal paths 16 . This for- 
malism is for instance used in classical diffusion driven 
systems where the transition rate from one stationary 
state to another, across a high barrier, is calculated. The 
path that crosses the barrier and has the largest probabil- 
ity is called the optimal path. Its shape and probability 
can be found using the principle of least action. Its prob- 
ability determines the transition rate. 

In this section, we will introduce the formalism of op- 
timal paths and apply it to our system, after which some 
generic properties of the action are derived. After this 
we show how the coherence time is related to the opti- 
mal paths. We investigate the dependence of the action 
of the optimal paths on the device parameters. This will 
be done both by using the results of previous section and 
by simulating the optimal paths. 



A. Trajectories and relation to Kramers' escape 
problem 

To start with, we outline the dynamics of the device 
in the absence of noise, describe the effect of noise and 
explain the relation to Kramers' escape problem. 

In Sec. |H] we have derived the requirements on stable 
stationary lasing in the HJL. Let us suppose that the op- 
tical field in the HJL is not in a point of stable lasing 
and that x — 0- If we disregard the noises we know that 
the field will evolve to one of the points. This evolution 
is described by Eq. ©. We visualize the dynamics by 
making a stream plot, as is done in Fig. [3l To describe 
the dynamics it is convenient to associate the field com- 
ponents, x, y, with coordinates of a 'particle' subject to 
a coordinate-dependent 'force field' that causes the mo- 
tion of the 'particle' with the velocity proportional to 
the force. Starting from a given initial condition, the 
'particle' will flow towards a stable stationary point. Its 
trajectory corresponds to a streamline in the plot. 
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FIG. 3. (Color Online) A stream plot corresponding to Eq. Q 
for variables x (horizontal) and y (vertical) , denned according 
to b = x + iy. The gray streamlines represent a 'force field'. 
When a 'particle' is placed in the 'force field' it will evolve 
along these lines to one of the attractors, indicated by A, de- 
pending on the initial condition. When placed at the thick 
red (dashed) or blue (not dashed) line, the 'particle' will not 
evolve to one of the attractors but to the saddle point, indi- 
cated by S. These lines form the separatrix of the system, 
that separates the domains of attraction of the attractors. 



Figure. [^corresponds to regime II as defined in Sec. ITT] 
In the figure the three stationary points are indicated by 
S, A and A'. The point S, at x, y = is the unsta- 
ble saddle point whereas the points A and A' are stable 
points or attractors. For a general initial condition in the 
vicinity of the saddle point the 'particle' is repelled from 
this vicinity. There is however a particular line called 
stable direction such that the 'particle' is attracted to 
the saddle point if the initial condition is chosen to be on 
this line. The streamlines coming to the saddle point at 
the stable direction form the separatrix of Eq. (j9)), which 
is the boundary separating the domains of attraction of 
the attractors A and A 1 . We see that the 'particle' cir- 
cles around an attractor before reaching it. This signifies 
that the 'force field' cannot be represented as a gradient 
of a scalar potential. The 'force field' must also have a 
transversal component, that causes the rotation and can 
be represented as the curl of a vector potential. 

Full equations of motion shall also include a variation 
of detuning \. With this we have three coupled equations 
of motion that define the 'force field' in the space of three 
variables 



x 



-x 



VA 2 - 1( X + x) - A + ti'i 



y, 



y = -y- 
x = -jx- 



VA 2 - 1( X + X) + A + Q"n 



(26) 



The last equation gives the relaxation of the detuning to 
its equilibrium position. The specifics of this equation is 
the presence of a saddle line The line x, y = is a saddle 
line because for each value of x the point x,y = is a 
saddle point with respect to the motion in x, y-directions. 

In Sec. IIIII noise was introduced. With the Langevin 
noise added to the dynamics, Fig. [3] no longer gives a full 



description of the 'particle' evolution. The 'particle' in 
the 'force field' will experience random kicks causing it 
to change from one trajectory to another, even when ini- 
tially at rest at an attractor. Effectively, the kicks cause 
the 'particle' to diffuse away from it's starting point, even 
against the streamlines of the 'force field'. Nonetheless, 
the 'force field' counters the diffusive flow trying to bring 
the 'particle' back to the attractor. These two kinds of 
motion have a fundamentally different nature. The evo- 
lution along the flow lines is deterministic, in contrast to 
the diffusive motion which is nondeterministic. Impor- 
tantly, the diffusive motion of the 'particle' lifts up the 
restriction of the 'particle' to stay in the domain of at- 
traction of an attractor. A series of random kicks may 
repel the 'particle' from an attractor and eventually pull 
it over the separatrix into the domain of attraction of the 
other attractor. After that the 'particle' relaxes naturally 
along a streamline to the other attractor. This trajectory 
corresponds to the large fluctuation we are after: a switch 
between the two stable lasing states has occurred. 

Let us suppose, in contrast to what is the case for 
Eq. ©, that a 'particle' moves in a potential well, of 
which the gradient defines the corresponding 'force field', 
and that the fluctuation-dissipation theorem is applica- 
ble, like in Kramers' escape problem. In this case the 
Fokker-Planck equation can be solved to calculate the es- 
cape probability from the potential wel l 15 ! 17 . This escape 
probability involves the Arrhenius factor, exp(W/T), 
where W is the height of the barrier that is crossed by 
the 'particle', counted from the potential minimum, and 
T is the effective temperature of the 'particle'. However, 
because of the transversal component of the 'force field' 
associated to Eq. © , we cannot use this general approach 
to solve Fokker-Planck equation (fTTj) . There, the escape 
probability does not depend on a potential barrier height, 
but on the shape and probability of the path that crosses 
the separatrix. We note that the potential barrier and 
the separatrix both form the boundary of the domain of 
attraction where the 'particle' is trapped. 



B. Optimal paths and the principle of least action 

With the conceptual framework in mind, let us con- 
tinue by showing how the probability to have a transition 
from one stationary state to the other, is related to the 
concept of the optimal path. We will only consider paths 
that start in stationary point A and end in stationary 
point A' . The time elapsing during the transition will be 
infinite, going from — oo to oo, because the 'particle' can 
only leave A and approach A' asymptotically. 

It is possible to express the probability distribution for 
the fluctuations in terms of an action. We start with a 
general set of Langevin equations, written down in vector 
form, 

?=V(F)+£(t). (27) 
The Langevin forces t;i(t) represent Gaussian white noise, 
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satisfying (&(t)) = and = A<M(* - t'), 

thus supposing an infinitesimally short autocorrelation 
time. The noise can be expressed as a single probability 
distribution that is a product of the probability distribu- 
tions of the separate Langevin forces at each moment in 
time, 



P = Cexp 



E 



2D, 



-dt' 



(28) 



with C the normalization constant. One may interpret 
this expression in another perspective, as the probabil- 
ity distribution for histories of fluctuations, because the 
integral in the exponent represents the history of fluctu- 
ations. Each history corresponds to a unique progress in 
time of £(£')• We stress that the values of £(i) at each 
time t remain strictly independent from the ones at other 
times t 7^ t. So how will this concept of histories aid us? 
Even though the Langevin forces might have infinitesi- 
mally short autocorrelation times, the vector field r will 
need some finite time to relax to the stationary state. 
During this time, the history of fluctuations will be rel- 
evant. The resulting cumulative effect of fluctuations is 
best captured by using Eq. (|27[) to rewrite the exponent 
of P, yielding 



Ce 



E- 



2D, 



dt' 



(29) 



Hence, we indeed see that the noise probability distribu- 
tion can be described using the action, defined as S. 

The probability distribution, P, is strongly reminis- 
cent of a path integral. Indeed, P is a probability density 
which, upon integration over all paths, gives the proba- 
bility to go from A to A'. When fluctuations are small, 
as is described in Sec. IIII1 the path integral is dominated 
by the 'classical' path and by the quadratic fluctuations 
about i1ji£. The 'classical' path is found by the principle 
of least action. It is, what we have called, the optimal 
path. The optimal path corresponds to the lowest value 
of the action for the transition, S op t, which relates to the 
transition probability 

P ~ F e ~ Sopt 

with F the prefactor corresponding to the quadratic fluc- 
tuations. Hence the probability of a transition from A to 
A' depends on the optimal path. 

Without specifying any details about V we note two 
things. First, by inspection we see that S > 0, with the 
equality satisfied for the deterministic path, for which 
f = V. Second, the integrand of S is equivalent to a 
Lagrangian and because it only implicitly depends on 
time, a conserved quantity similar to the Hamiltonian 
exists 



= 0. 



(30) 



The equality of the Hamiltonian to zero is not univer- 
sal, but is appropriate for the path we consider. This 
path connects the stationary points, where r = V = 
and therefore H — 0. Being a conserved quantity, H 
must always be zero at the path connecting the station- 
ary points. 

The action as it is given in Eq. ([2TR) is not in the most 
convenient form. To minimize this action we need to find 
both the shape of the optimal path and the speed with 
which the path is travelled. By exploiting the conserved 
Hamiltonian, we can relate the speed to the 'force field' 
resulting in an action where the optimal path is only de- 
termined by its shape. To show this, we choose the diffu- 
sion constants, without loosing generality, as Di = 1 for 
all i. This is allowed, because we can always rescale 
and Vi accordingly. Applying this to Eq. (|30|) we then 
find that |r| = \V\. With this, the terms r 2 /2 and 
V 2 /2 can be rewritten and added to each other to be- 
come |r||V|. Because both terms in the Lagrangian now 
only contain a single time derivative, we can parameterize 
time at will. We parameterize it as a monotonously grow- 
ing function ('(t) > 0, such that C(*o) = and = 1. 
Consequently f, = ((df/dQ, so that the action becomes 



S = 



Of 



\V\ -^-V\<l< (■'!!! 



Minimizing this action will result in the optimal path 
irrespective of the exact form of parametrization £(t) and 
therefore irrespective of the speed with which the optimal 
path is travelled. 

From the expressions for the action, Eq. (|31[) . an in- 
tuitive understanding of the dynamics of the system can 
be obtained. Suppose that we are at some point in the 
'force field' r and move away from this point a random 
infinitesimal distance dr. The probability to move in a 
certain direction is then determined by the increase in 
action associated to this direction, which is in turn de- 
termined by the angle between the 'force field' vector, V 
and dr. With the two vectors aligned the action is small- 
est (zero) and the probability to move in that direction 
the highest. Again, this corresponds to the determinis- 
tic solution. With the vectors counter aligned the action 
is largest and the probability to move in that direction 
the lowest. Close to stable or unstable stationary points, 
when \r\ 2 is of the order of the diffusion constants, the 
absolute difference between alignment and counter align- 
ment is small because \V\ is close to zero. The probabil- 
ities to go in any direction become of the same order of 
magnitude. This is the regime where the quadratic fluc- 
tuations are important. Far away from stationary points, 
however, | V| can be quite large so that the absolute differ- 
ence between alignment and counter alignment becomes 
large. Then, because the action appears in the exponent, 
the probability to follow the deterministic path will be 
very close to unity. Going in some other direction will 
cause a large increase in action. 

Concluding this subsection we make two remarks. 
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First, the probability distribution of Eq. ((29)) is a solu- 
tion to the Fokker-Planck equation (fTTj) in the Wentzel- 
Kramers-Brillouin approximatio n 18 ! 19 . The solution is 
exact in the limit where the diffusion constants are in- 
finitesimally small. Second, the coherence time, or escape 
time, is related to the transition probability^, yielding 



Ke- 



rn 



Here, K relates to the quadratic fluctuations around the 
optimal path. We will disregard this prefactor and only 
consider the action of the optimal path, which determines 
the order of magnitude of the coherence time. 



C. The action for the HJL 

Let us apply the above described ideas and calculate 
the action for the HJL, using the Fokker-Planck equation 
of Eq. pT|) . The Fokker-Planck equation contains three 
variables. Two of them, x,y, associated to the optical 
field, defined as b = x + iy, and one associated to the 
fluctuations in detuning, \. We redefine the latter as z = 
x/ \JDy_-, where = 27 ((x 2 )) is the diffusion constant 
of x- As a result, the diffusion constants associated to 
fluctuations in x, y, z are all equal to one. Because any 
Fokker-Planck equation, with constant diffusion terms, 
has an equivalent Langevin equation 15 , Eq. (Ilip can be 
written in the form of Eq. (|27]l , with 



(-x + x V^~T y + ^ 

V= \-y- X ^W^\x-\f- 

U = ft + y/2-y((x 2 ))(A 2 - 1) z(x 2 + y 2 ) (33) 

+ io'V + y 2 ) 2 + ^ 2 -2/ 2 ) 

D x = D y = D z = l 

where f= (x,y,z) a dimensionless quantity. We remind 
the reader that A, 7, w, Q,, SI', SI" and U [Eq. ©] are 
dimensionless with the corresponding frequencies being 
measured in units of T/2. The drift terms V can be 
explicitly written in terms of scalar and vector potentials 



V = + V x A, 



/ 2 1 2\ , 7 2 

oO +y ) + 77 z 



(34) 



A = A z z, A z = -[u>(x 2 + y 2 ) + U], 



demonstrating that the vector field corresponding to V 
contains both longitudinal and transversal parts. The 
action for the HJL is now simply given by Eq. (j3TJ) with 
the vector and the diffusion constants given by Eq. 



D. Estimation of the action 

Now that we have discussed the formalism of optimal 
paths and have established the relation with the escape 



time, we want to find explicitly the value of the optimal 
action, depending on the system parameters of the HJL. 
Before that, let us present the estimation of the action 
for the optimal path, S op t, in different regimes. For this 
we will use the result of Sec. IIIII In the following subsec- 
tion we will simulate numerically the optimal paths and 
the corresponding action, to validate the findings of this 
subsection. 

In Sec. |HT] we have discussed the fluctuations in the 
device in the linear noise approximation. Under these 
assumptions the fluctuations are small and obey a Gaus- 
sian distribution with variances given by Eq. (1201) . 

For the order-of-value estimations we can use this for- 
mula for small fluctuations to find the probability of large 
fluctuations. In particular, large fluctuations in n are im- 
portant, because at n = 0, the 'particle' reaches the sad- 
dle point enabling the switching to the opposite stable 
point. Large fluctuations in either of the other variables 
will not drive the 'particle' across the separatrix. The 
probability for this fluctuation in Gaussian approxima- 
tion is simply proportional to exp[— n 2 /(2(<5n 2 ))], with 
(Sn 2 ) given by Eq. (|2"0|) 



(Si 



2( 7 + 2)(A 2 -l)(x 2 ) 



(l + X )[4(-A 2 -l)(l + x)+7 2 + 2 7 ] 
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2n, 



1 



(Ai-l)(l + x) 



(35) 



We can thus estimate S op t — n 2 /(8n 2 ), where the co- 
efficient is different in different parameter regimes. Let 
us list here all possible limiting cases corresponding to 
different relations between parameters of the device. 

In the case (i) the large fluctuation is caused by quan- 
tum photon noise. The fluctuation is estimated by the 
last term in the expression (|35j) and the optimal action 
is given by 



S opt = n s (l+x)Fi(A, x) 



(36) 



F\ being a coefficient of the order of one. The probabil- 
ity of the large fluctuation is enhanced near the lasing 
threshold, x = — 1, and is exponentially small otherwise, 
provided n s ^> 1. 

In all other cases the probability of the large fluctu- 
ation is mainly caused by the detuning noise. In the 
limit 7 <C 1, when the switching time of the emitters is 
longer than the damping time of the optical resonator, 
we can distinguish two cases. The case (ii) takes place 
near the lasing threshold, Awl, and is characterized by 
7 3> A 2 — 1. The estimation for S opt reads 



7(1 + x) 



i)(r 



-Fiix) 



(37) 



and is inversely proportional to the variance of the de- 
tuning fluctuations. In the opposite case (iii) 7 -C A 2 — 1. 
The estimation is given by 



(1 + x) 2 

Sopt — ^~2^ ^3 (AX) 



(38) 



16 



and does not depend on 7. 

If the switching times are relatively short, 7 > 1, we 
distinguish two cases depending on the ratio of 7 and A. 
In case (iv) 7 3> A and the estimation reads 

S° P t = ^y^(x) (39) 
while in the opposite case (v) (7 <C A) it reads 

Sopt = { -^-F 5 ( X ). (40) 

7 Or) 

In all these cases coefficients F2-5 are supposed to be 
of the order of one. The boundaries in parameter space 
separating the case (i) from the cases (ii)-(v) are found 
by comparing the corresponding estimations. 



E. Numerical results 

Now that we have the estimations of the optimal ac- 
tion and the understanding of its dependence on relevant 
parameters, we quantify numerically the shape of the op- 
timal path and the action. We describe the numerical 
method and we comment briefly on its limitations. We 
compare the numerical results and the estimations and 
discuss the peculiarities of the path shapes. 

We use the action of Eq. (f^Tj) . A path f (£) starts at the 
attractor at positive x, negative y and z = (Am Fig. [3]) 
and ends at the point A'. The path can be separated in 
two parts. During the first part the path proceeds from 
A till a point at the separatrix that divides the domains 
of attraction. During the second part the path follows 
the natural motion of the 'particle' satisfying equations 
of motion Eq. (f2"6")l . This part does not contribute to the 
action, since dr/dC, \\ V, and thus does not have to be 
evaluated numerically. In all cases investigated numeri- 
cally, we have found that the first part of the path end 
at the saddle line x, y — at some yet to be optimized 
value of z. 

To find the optimal path we could first formulate a 
boundary value problem, using the differential Euler- 
Lagrange equation from the action Eq. (l3"3")l . The solution 
to this is the optimal path, which can only be found using 
a shooting method. With this method, one changes ini- 
tial conditions and integrates the Euler-Lagrange equa- 
tion trying to reach the predefined final condition. This 
method is not efficient to find the trajectories we are look- 
ing for because they have a relatively low probability. 

We implement a more efficient method. We directly 
minimize the action of a predefined path that depends 
on a finite number of parameters as function of the path 
parameters. For the path parametrization, we implement 
a Bezier curve of some sufficiently high order n (from 12 
to 16 depending on the region in the parameter space) 



(41) 



The order plus one equals to the number of points Pi in 
xyz space that define the curve. The starting point of the 
Bezier curve, Po, is fixed to the attractor while the last 
point, P n is fixed to the saddle line, P* being the subject 
of the optimization procedure, as well as the coordinates 
of all other points. The accuracy of this method depends 
on the order of the Bezier curve used. The higher the or- 
der, the better the accuracy. The order required to have a 
sufficiently accurate path will depend on the complexity 
of the shape of the path. 

We have calculated the optimal paths for several pa- 
rameter regimes, to validate the estimations of Eqs. (|36|) - 
(|4"0"|) and find numerical values of the prefactors. In 
Fig. BJa) we show the dependence of the value of the op- 
timal action on the number of photons in the resonator 
mode, n, for three different values of 7 and x — 0- F° r 
sufficiently small values of n, S is linear in n and the 
case (i) applies, yielding F% w 0.4 for all three values of 
7. For large values of n the action S saturates to a value 
depending on 7. For 7 = 0.1 the case (iii) applies, with 
F3 « 0.5, while the case (iv) applies for 7 = 1000 with 
F4 1=3 0.5. The curve for 7 = 10 at large n corresponds to 
a crossover regime between the cases (iv) and (v). 

In Fig.HJb) we show the dependence of the value of the 
action on 7 for two values of n(x 2 )- In the limit 7^1 
case (iii) applies, with also here F3 ~ 0.5. In the limit 
7 » A case (iv) applies with F 4 depending on (x 2 )- At 
n(x 2 } = 25 F4 fa 0.5, which is consistent with its value 
found previously, whereas at n(x 2 ) = 100 F 4 rj 0.8. 

Figure [He) shows the dependence of the value of the 
action on the detuning x f° r three values of n(x 2 )- For 
n(x 2 ) = 0.25 case (i) applies with Fi ~ 0.4, which is 
consistent with the value of F\ for Fig.^Ja). For n(x 2 ) — 
25 case (iii) applies with F 3 ss 0.5, which is also consistent 
with the value previously found. The curve with n(x 2 ) = 
1.5 corresponds to the crossover regime between cases (i) 
and (iii) so that no value can be extracted for any of 

We study the shapes of the optimal paths for the 
cases listed, with representative examples given in Fig. [5] 
Rather surprisingly, we find that the shapes of the op- 
timal paths are qualitatively similar to the shapes of 
the deterministic paths, that are defined by Eqs. (|26l) . 
Both the optimal path and the deterministic path circle 
counter clockwise around the stable point / while mov- 
ing respectively from/to the point. The three optimal 
paths in the right column of Fig. [5j with n(x 2 ) = 0.03, 
correspond to the case (i). In this case, the intrinsic 
photon number fluctuations dominate. As a signature 
of this, the path circles around the stationary point / in 
the xy-plane, without any significant displacement in the 
z-direction. 

The optimal paths in the column with n(x 2 } = 30 cor- 
respond to the regime where the fluctuations in detuning 
dominate. Here the displacements in the z-direction are 
two orders of magnitude larger than for the optimal paths 
of the previous case. The upper left path corresponds to 
the case (iii), where the photon number fluctuations are 
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FIG. 4. (Color Online) The value of the action of the optimal path for several parameter regimes. Fixed parameters are denoted 
in the graphs. Parameters were chosen such that we are far away from the lasing threshold where the use of optimal paths is 
justified. The main text contains a comparison with Eqs. (|36|) - (|40|) to estimate the values F1-5. (a) The action as a function 
of n(x 2 ), for three values of 7. The axes are logarithmically scaled. For small values of n(x 2 ) the action is linear in photon 
number and independent of 7. This is the regime where intrinsic photon number fluctuations dominate. For large n(x 2 ) the 
detuning fluctuations dominate and the action saturates to a value depending on 7. (b) The action as a function of 7 for two 
values of n(x 2 }, with logarithmically scaled axes. For 7 <?C 1 the (x 2 )S is constant while for 7 3> A it is proportional to 7. In 
the region between we see a minimum at 7 ~ A. The curve for n(x 2 ) = 25 has points in common with the three curves of (a), 
(c) The action as a function of the detuning x, for three values of n(x 2 }- 



thus fast that the optical field adiabatically follows the 
fluctuations in the detuning. The path has a shape re- 
sembling a piece of parabola and reaches the saddle line 
close to the value of z corresponding to the lasing thresh- 
old (Fig. [TJ. For the path with 7 = 1000 [case (iv)], de- 
tuning fluctuations occur on shorter timescales so that 
the photon number fluctuations become relatively more 
important, leading to an increased number of circulations 
of the path around the stationary point I. The other op- 
timal paths correspond to crossovers between the cases. 
The general trend is however clear: with increasing 7 
(n(x 2 )) the displacement of the optimal path in the z- 
direction (xy-plane) becomes smaller. 

We briefly comment on the accuracy of the method 
used. The order of the Bezier curve that is required to 
have an accurate optimal path is different for the different 
cases. The simple optimal path of case (iii) is found with 
high accuracy (up to one percent) with a Bezier curve of 
the order 12, while the optimal path with 7 = 1000 and 
n(% 2 ) = 1.5 is a Bezier curve of order 14 and is less ac- 
curate (up to ten percent). We do not use a higher order 
because it would greatly increase the computation time 
without providing new insight. The shape of the opti- 
mal path can however be inferred, especially for the case 
(i). Since the optimal path resembles the deterministic 
path and the deterministic path in the case (i) circles 
~ A/2tt around the stable point (Sec|TT]), we expect that 
also the optimal path circles ~ A/2ir times around I. To 
calculate this optimal path requires a Bezier curve of at 
least the order of A/2ir. The method is thus inefficient 
for very big A. For the value A = 20 3> 1 that was 
used, we are far away from the lasing threshold, while 
being able to find the optimal path with sufficient accu- 
racy. For the cases other than the case (i) the number 
of times the path circles around / decreases, resulting in 
more accurate optimal paths. 



V. CONCLUSIONS 

In this Article we have shown that a device that com- 
bines stimulated emission of light with superconductivity, 
can produce coherent radiation with exponentially long 
coherence times. 

We have introduced a general model for such Half- 
Josephson Laser (HJL), where a big number of quan- 
tum states are coupled to two superconducting leads and 
are capable of emitting light by electron-hole recombi- 
nation into a resonant optical mode. We argued that 
there should also be incoherent emission of photons in 
the environment that creates an imbalance in the popula- 
tion of emitter states required for lasing. The stochastic 
dynamics of the HJL is compactly described by a sin- 
gle Fokker-Planck equation, Eq. (fTTj) , for three real vari- 
ables, incorporating two noise sources. One noise source 
is the intrinsic quantum noise associated to the discrete 
changes of photon number in the resonator mode. The 
other noise source is related to the random switching be- 
tween the eigenstates of the quantum emitters and causes 
the fluctuations of detuning. 

We have found a lasing regime where the nonlasing 
state is unstable. The optical phase is locked to the su- 
perconducting phase difference in two possible stable las- 
ing states with a phase difference of tt. In this case, small 
fluctuations do not lead to decoherence, which is caused 
by a large fluctuation whereby the device switches be- 
tween the two stable lasing states. To evaluate the coher- 
ence time of the HJL, we study these large fluctuations 
estimating their probability by the optimal path method. 
Our numerical results show agreement with simple esti- 
mations deduced from the distribution of small fluctua- 
tions. 

We have found the decoherence times to be exponen- 
tially long ~ exp[— 5 op t] provided 5 op t ^ 1- The optimal 
action is inversely proportional either to the number of 
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FIG. 5. Optimal paths corresponding to points in Fig. |4]for various values of n(x 2 } and 7 as indicated above the columns and 
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photons in the laser or to the variance of the detuning. 
The latter is inversely proportional to the number N of 
emitters in the laser. Thereby we prove the feasibility of 
the exponentially long times for a big number of emitters. 
S op t was a ^ so sri0wn to depend on 7, the ratio of emitter 
switching rates to damping rate from the resonant mode, 
and the distance to the lasing threshold. 

To give an example we can express these conditions 
for a simple microscopical model formulated in terms of 
the parameters of the interaction Hamiltonians [Eqs. ^ 
and fll}] hi Sec. [I] Firstly, we estimate the system pa- 
rameters A and fl" [Eq.flSJ] assuming that the energy of 
the states of the emitters is only perturbed by these two 
interactions. Taking E s as the typical energy scale of the 
eigenstate of a single emitter and using perturbation the- 
ory we can therefore estimate the order of magnitude of 
A and O" 



A~ N 



\A e A h \G 2 



n" 



A- 



(42) 



Here it is assumed that the population of emitter states 
is asymmetric as discussed in Sec HI To have lasing in the 
HJL requires A > T. Additionally, a large number of pho- 
tons, n, in the resonator mode requires A 3> f2" implying 
that I A e A/j I 3> G 2 . Secondly, we estimate the variance 
of the detuning fluctuations, x, which correspond to the 
fluctuations of f2' [Eq. ©] scaled with A. Using pertur- 
bation theory, we estimate (f2') ~ NG 2 /E S . Its variance 
is of the order of a(f2') 2 /N , with a < 1 being a constant 
depending on the transition rates between the emitter 



eigenstates, where it is of the order of one when all rates 
are of the same order. This leads to 



(x 2 



Et 



~N\A e A h 



(43) 



where the estimation for A was used. Typically, this 
quantity should be much smaller than one. 

We list several ideas for possible applications of the 
HJL. In the HJL model, we neglected voltage fluctua- 
tions of the bias. When present, these fluctuations will 
transfer to fluctuations in the optical phase of the laser. 
By measuring these and providing negative feedback the 
voltage fluctuations, and hence the optical fluctuations, 
can be reduced and stabilized. Another idea is similar to 
what was proposed in RefH One can embed two HJLs 
in a SQUID and let them interact optically, so that an 
Aharonov-Bohm effect manifests in a shift of the optical 
phases that is flux-dependent. As an alternative to the 
latter idea, one can also let two distant SQUIDs interact 
via optically coupled HJLs. This may allow for coherent 
interactions between widely separated superconducting 
qubits. Finally, the exponentially long coherence time 
of the HJL can make it very suitable for metrology pur- 
poses, where the measurement of optical frequencies in 
the next-generation atomic clocks require ultrastable op- 
tical lasers^. 
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